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We calculate the masses of bottom mesons using an improved relativistic action for the 6-quarks 
and the RBC/UKQCD Iwasaki gauge configurations with 2+1 flavors of dynamical domain-wall 
light quarks. We analyze configurations with two lattice spacings: a~ = 1.729 GeV (a w 0.11 fm) 
and a -1 = 2.281 GeV (a « 0.086 fm). We use an anisotropic, clover-improved Wilson action for the 
6-quark, and tune the three parameters of the action nonperturbatively such that they reproduce the 
experimental values of the B s and B* heavy-light meson states. The masses and mass-splittings of 
the low-lying bottomonium states (such as the r\b and T) can then be computed with no additional 
inputs, and comparison between these predictions and experiment provides a test of the validity 
of our method. We obtain bottomonium masses with total uncertainties of ~ 0.5-0.6% and fine- 
structure splittings with uncertainties of ~ 35-45%; for all cases we find good agreement with 
experiment. The parameters of the relativistic heavy-quark action tuned for fe-quarks presented 
in this work can be used for precise calculations of weak matrix elements such as B-meson decay 
constants and mixing parameters with lattice discretization errors that are of the same size as in 
light pseudoscalar meson quantities. This general method can also be used for charmed meson 
masses and matrix elements if the parameters of the heavy-quark action are appropriately tuned. 



PACS numbers: 

I. INTRODUCTION 

Precise knowledge of mass spectrum, decay, and mix- 
ing properties of hadrons containing one or more bottom 
or charm quarks is essential to advancing our understand- 
ing of the parameters of the Standard Model. Lattice 
Quantum Chromodynamics (QCD) provides methods to 
compute these quantities from first principles. Conven- 
tional lattice calculations with heavy quarks are challeng- 
ing, however, because it is impractical to use a sufficiently 
small lattice spacing to control the 0(ma) n discretization 
errors directly. 

One way to address this challenge is to adapt the lat- 
tice description of heavy quarks to correctly describe 
heavy-quark physics in a carefully circumscribed kine- 
matic range. Such approaches include heavy-quark ef- 
fective theory (HQET) [1] in which the limit of infinite 
quark mass is considered and the continuum limit of the 
lattice theory reproduces the continuum heavy quark ef- 
fective theory. A second method is non-relativistic QCD 
(NRQCD) [2, 3] in which the mass of the heavy quark is 
assumed to be much greater than the inverse lattice spac- 
ing but the momentum-dependence of the heavy quark 
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energy is included in the non-relativistic limit. Each of 
these approaches has its own limitations. Specifically, ra- 
diative corrections to the coefficients of the NRQCD La- 
grangian contain power-law divergences that blow up in 
the limit ma —> 0, while HQET cannot deal with quarko- 
nia. 

The Fermilab or relativistic heavy quark (RHQ) 
method [4-6] provides a more complete framework for 
heavy-quark physics. It applies for all values of the heavy 
quark mass tuq, for both heavy- heavy and heavy-light 
systems, and allows a continuum limit. The improved 
RHQ action accurately describes energies and amplitudes 
of on-shell states containing heavy quarks whose spatial 
momentum p is small compared to the lattice spacing. 
It can be shown [6] that all errors of order \pa\, (toqo)™ 
and \pa\{mQa) n for all non-negative integers n can be re- 
moved if an anisotropic, clover-improved Wilson action 
is used for the heavy quark. This action depends on 
three relevant parameters: the bare quark mass too, an 
anisotropy parameter ((mod) and the coefficient cp(moa) 
of an isotropic Sheikholeslami and Wohlert term. 

In order to exploit this RHQ approach, values for these 
three parameters are needed. The bare charm or bot- 
tom quark mass, to , must be determined from exper- 
iment, usually by equating the known mass of a phys- 
ical state containing one or two heavy quarks with the 
mass determined from a lattice calculation with the RHQ 
action. The remaining two parameters, £ and cp, may 
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be estimated from lattice perturbation theory or deter- 
mined with a nonperturbative technique. We cannot use 
the nonperturbative method of Rome-Southampton [7] to 
tune cp and £ because the Rome-Southampton approach 
depends on evaluating off-shell amplitudes, whereas the 
3-parameter RHQ action only controls discretization er- 
rors for on-shell states. On-shell step-scaling methods 
can be used, either via the Schrodinger functional ap- 
proach or a simple comparison of small volume spectra 
between calculations with varying lattice scale but iden- 
tical physical volumes [8] . Both of these step-scaling ap- 
proaches, however, involve substantial computational ef- 
fort, requiring a scries of carefully matched finite volume 
simulations with varying lattice spacing. 

In the work reported here, we use another approach 
and determine the two remaining parameters £ and cp 
nonperturbatively by imposing two simple conditions. 
The first condition is the often-exploited requirement 
that the energy of a specific heavy-heavy or heavy-light 
state depend on that state's spatial momentum in a 
fashion consistent with continuum relativity: E(p) 2 = 
p 2 + A/ 2 . The second constraint is that a specific mass 
splitting agree with its experimental value. For the case 
of bottom, a natural choice is the B* — B s mass splitting. 
Thus, using the bottom system as an example and in- 
cluding the bare quark mass mo, we determine our three 
parameters mo, £( m oa) and cp(moa) by requiring that 
experimental values are obtained for mp s and mp* and 
that Ep s has the proper dependence on pp s . 

As is described below, these three conditions are 
straightforward to impose and yield quite precise results 
for the three unknown parameters. This approach has 
the disadvantage that a possible experimental predic- 
tion from lattice QCD, a non-trivial spin-spin splitting, 
cannot be made. With this approach, however, we can 
immediately determine the masses of a large number of 
heavy-heavy and heavy-light states. These results can 
be viewed as tests of QCD and can be used to explore 
the accuracy and limitations of the RHQ approach. Fi- 
nally, once the RHQ action has been determined by fix- 
ing these three parameters, it can be used to compute 
phenomenologically-important charm and bottom decay 
constants and mixing matrix elements, which are needed 
for determinations of CKM matrix elements and con- 
straints on the CKM unitarity triangle 

In this paper we present results for the bottom sys- 
tem. Our calculation is performed on the 2+1 flavor, 
domain wall fermion (DWF) + Iwasaki gauge-field en- 
sembles generated by the LHP, RBC, and UKQCD col- 
laborations with several values of the light dynamical 
quark mass at two lattice spacings, a ps 0.11 fm and 
a Rj 0.086 fm [9, 10]. For the heavy-light mesons, the 
heavy quark will typically carry a small spatial momen- 
tum, \p\ pa Aqcd- Thus, for these systems the expected 
|pa| 2 errors are of the same size as those encountered in 
calculations involving only light quarks. For heavy-heavy 
systems, however, the spatial momenta will be larger: 



\p\ ps a s rriQ, where ttiq is the heavy-quark mass and a s 
the strong interaction coupling constant evaluated at a 
scale appropriate for such a bound state. While for char- 
monium a s mQ may be on the order of Aqcd, this is not 
the case for bottomonium where discretization errors are 
expected to be three to four times larger due to the larger 
b-quark mass (mt,/m c ps 3.3). Thus we choose to tune 
the RHQ action for 6-quarks using bottom-light states in 
order to minimize systematic uncertainties. In particu- 
lar, we match to the experimentally-measured masses of 
the bottom-strange states B s and B* in order to avoid 
the need to perform a chiral extrapolation in the valence 
light-quark mass. Once we have determined the values 
of the parameters {mo,cp,£} we make predictions for 
the masses and mass-splittings of several low-lying bot- 
tomonium states: t]b, T, \bo, Xbi, and hb- Our results 
agree with experiment within estimates of systematic un- 
certainties, confirming the validity of the RHQ approach 
and bolstering confidence in future computations of weak 
matrix elements with the RHQ action. 

This work was begun by Li, who presented prelimi- 
nary values for the RHQ parameters and bottomonium 
masses on the coarser a ps 0.11 fm ensemble at Lattice 
2008 [11]. We improve upon his results in several ways, 
most notably in determining the RHQ parameters solely 
from quantities in the bottom-strange system. (Li used a 
hybrid of bottom-strange and bottomonium observables 
for the tuning.) This reduces the systematic errors in 
the resulting parameters due to heavy-quark discretiza- 
tion effects, as discussed above. We also significantly 
increased the statistics, more than quadrupling the num- 
ber of configurations analyzed, and optimized the spa- 
tial smearing wavefunction used for the 6-quarks in or- 
der to reduce excited-state contamination in the bottom- 
strange 2-point correlators. More recently Peng extended 
this work to the finer a ps 0.086 fm ensembles and pre- 
sented preliminary values for the RHQ parameters and 
bottomonium masses at Lattice 2010 [12]. Again, we 
polish this result with increased statistics and improved 
b-quark smcarings. 

This paper is organized as follows. In Section II we first 
present the form of the relativistic heavy-quark action 
used in this work. We then describe the numerical strat- 
egy used to determine the three parameters mo, £(moa) 
and cp(moa). Next, in Section III we present the tuning 
of the RHQ parameters for bottom. We give the actions 
and parameters used in our numerical simulations, and 
then discuss the fits of heavy-light meson 2-point correla- 
tors needed to extract the lattice values of the B s and B* 
meson masses. Using this data we tune the parameters of 
the RHQ action such that it applies to 6-quarks. In Sec- 
tion IV we use the resulting RHQ parameters to predict 
the masses of several bottomonium states and compare 
the results with experiment. Finally, we summarize our 
results and discuss future plans in Section V. 
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II. FRAMEWORK OF THE CALCULATION 

A. Heavy-quark action 

The relativistic heavy-quark method provides a consis- 
tent framework for describing both light quarks (arriQ <C 
1) and heavy quarks (amo <; 1) [4-6]. This approach re- 
lies upon the fact that, in the rest frame of bound states 
containing one or more heavy quarks, the spatial momen- 
tum carried by each heavy quark is smaller than the mass 
of the heavy quark: for heavy-heavy systems |p| ~ a s mo 



and for heavy-light systems \p\ ~ Aqcd- Then one can 
perform the usual Symanzik expansion in powers of the 
spatial derivative -D; (which brings down powers of ap). 
One must, however, include terms of all orders in the 
mass mod and the temporal derivative Dq. This sug- 
gests that a suitable lattice formulation for heavy quarks 
should break the axis-interchange symmetry between the 
spatial and temporal directions. 

In this work we use the following anisotropic clover- 
improved Wilson action for the 6-quarks: 



where 
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and 7/i = 7 t , {7^,7,,} = 25^ and a^ v = §[7^,7,,]. 
Christ, Li, and Lin showed in Ref. [6] that one can absorb 
all positive powers of the temporal derivative by allow- 
ing the coefficients cp(m§a) and ((mod) to be functions 
of the bare-quark mass niga. Thus, by suitably tuning 
the three coefficients in the action - the bare-quark mass 
mo a, anisotropy parameter £, and clover coefficient cp 
- one can eliminate errors of 0(|p|et), O([moa]"), and 
O(\pa\[moa] n ) from on-shell Green's functions. The re- 
sulting action can be used to describe heavy quarks with 
mod, > 1 with discretization errors that are comparable 
to those for light-quark systems. 

The relativistic heavy quark formulation based on 
Ref. [6] and used in this work is one of several vari- 
ations. This general approach was first introduced by 
El-Khadra, Kronfeld, and Mackenzie in Ref. [4], and has 
been used recently by the Fermilab Lattice and MILC col- 
laborations for many phenomenolgical applications such 
as decay constant and form- factor computations [13-16]. 
In practice, however, Fermilab/MILC use a different ap- 
proach to tune the parameters in the action, Eq. (1), 
than our method described below in Sec. II B. They fix 
the anisotropy parameter £ to unity and the clover coeffi- 
cient cp to its tree-level mean-field improved lattice per- 
turbation theory value (1/wq), and then tune only the 
hopping parameter k (which is equivalent to the bare- 
quark mass) nonperturbatively [17]. The Tsukuba group 
uses a slightly different formulation of the action in which 



they do not use field rotations to eliminate redundant 
operators [18]; hence their version of the action has four 
unknown coefficients rather than three in the RHQ or 
Fermilab variants. For on-shell Green's functions the 
Tsukuba and RHQ/Fermilab actions are equivalent. In 
practice, however, the inclusion of redundant couplings 
means that one cannot nonperturbatively tune all four 
parameters simultaneously by only adjusting the energies 
of heavy hadrons because one will run into flat directions 
in parameter space, as was shown in Ref. [6] . Hence they 
rely upon lattice perturbation theory for quark-quark 
scattering amplitudes to determine at least one of the 
coefficients [18]. 

Because the lattice action breaks Lorentz symmetry, 
mesons receive corrections to their energy-momentum 
dispersion relation due to lattice artifacts: 

{aEf = (aAhf + (ap) 2 + 0{[apf) . (5) 

The quantities Mi and M2 are known as the rest mass 
and kinetic mass, respectively, 

/ dE 2 \ 1 

M 1 =E(p = 0), M 2 =MtX (— T ) , (6) 

and will generally be different for generic values of the 
parameters {moa,cp,£}. We will exploit this fact in the 
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RHQ tuning procedure described in the following subsec- 
tion. 



B. Parameter tuning methodology 

We tune the values for the RHQ parameters 
{moa, cp, to describe bottom or charm quarks by re- 
quiring that calculations of specified physical on-shell 
quantities with the action in Eq. (1) correctly reproduce 
the experimentally-measured results. In particular, for b- 
quarks we determine the RHQ action using the bottom- 
strange system because both discretization errors and 
chiral extrapolation errors are expected to be small. We 
match to the experimental values of the spin-averaged B s 
meson mass: 



M Bs = X - {M Bs 



3M B .) 



and hyperfine splitting: 



AMr = M F 



(7) 



(8) 



We also require that the B s meson rest and kinetic masses 
are equal: 



Mt 



Mk 



(9) 



so that the B s meson satisfies the continuum energy- 
momentum dispersion relation E B (p) = Pb + M% ■ 
(Note that throughout this work we denote meson masses 
with a capital "M" and quark masses with a lower-case 
"to" in order to avoid confusion in situations where the 
context is insufficient.) We could in principle have used 
other states (e.g. scalar or vector mesons), other mass- 
splittings (e.g. the spin-orbit splitting), or even other sys- 
tems (e.g. heavy- heavy mesons) to tune the parameters 
of the RHQ action, since the same parameters should de- 
scribe 6-quarks in all of these arenas. Instead, however, 
we can make predictions for these quantities using the 
tuned values of {moa,cp,£} and use them to test the 
validity of our approach. 

We determine the tuned values of {tooo, cp, £} nonper- 
turbatively using an iterative procedure. The bottom- 
strange meson masses in general will have a nonlinear 
dependence on the RHQ parameters. We choose to work 
in a region sufficiently close to the true parameters such 
that the following linear approximation holds: 
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where J is a 3 x 3 matrix containing the linear coef- 
ficients (analogous to the slope in the lxl case) and 
A is a 3-element column vector containing the con- 
stants (analogous to the intercept). For a single step 



of the iteration procedure we compute the quantities 
{Mb, , AMfl a , M-^ s /M B "} at seven sets of parameters 
(see Fig. 1) in which we vary one of the three parameters 
{mod, cp,C} by a chosen uncertainty ±<j{ mga cp q (not 
to be confused with the statistical errors in the tuned pa- 
rameters {tooo, cp, £}) while holding the other two fixed: 
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This allows us to test whether or not the "box" of param- 
eter space defined by the seven parameter sets in Fig. 1 is 
in the linear region such that Eq. (10) applies. If indeed 
we are in the linear region, we then compute the matrix J 
and vector A via a simple finite difference approximation 
of the derivatives: 



J = 



Y 3 -Y 2 Y 5 ~Y 4 Y 7 -Y 6 



2cr„ 
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A = Y\ — J x [mod, cp, C] 



(12) 
(13) 



where Yi is the 3-element column vector containing the 
values of meson masses and splittings measured on the 



;th 



parameter set listed in Eq. (11) 
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Finally, the tuned RHQ parameters are given by: 
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(15) 



We consider the RHQ parameters {TOoa, cp,C} to be 
"tuned" when all three of the values obtained via Eq. (15) 
are within the "box" defined by the seven parameter sets 
in Fig. 1. This condition ensures that we are interpo- 
lating, rather than extrapolating, to the tuned point. If 
the result for any of the parameters {mo a, cp, £} lies out- 
side the box, we re-center the box around the result of 
Eq. (15) and perform another iteration step. We repeat 
this procedure until all three tuned RHQ parameters lie 
inside the box. 

Once the RHQ parameters have been tuned, we can 
use them to predict the masses of other heavy-light and 
heavy-heavy meson states, and ultimately to compute 
heavy-light meson weak matrix elements. We compute 
the desired quantities on the same seven sets of parame- 
ters used for the final iteration of the tuning procedure. 
We then propagate the statistical errors in the tuned 
RHQ parameters to these quantities using the jackknife 
method; this accounts for correlations between the pa- 
rameters toq, cp, and £. 



5 




3.5 
3.3 
3.1 
2.9 



8 ^°8-30li5 sTe c 5 - 8 

moa cp 

FIG. 1: Location of the seven sets of parameters used to 
obtain the tuned values of {moa, Cp, £}• 



III. LATTICE CALCULATION OF RHQ 
PARAMETERS FOR BOTTOM 

A. Lattice simulation parameters 

The parameters of the RHQ action suitable for de- 
scribing 6-quarks depend upon the choice of actions for 
the gauge fields and sea quarks. In this work we perform 
our numerical lattice computations on the "2+1" flavor 
domain-wall fermion ensembles generated by the LHP, 
RBC, and UKQCD Collaborations [9, 10]. These lat- 
tices include the effects of three light dynamical quarks; 
the lighter two sea quarks are degenerate and we denote 
their mass by to;, while the heavier sea quark, whose 
mass we denote by to/,, is a little heavier than the phys- 
ical strange quark. The RBC/UKQCD lattices com- 
bine the Iwasaki action for the gluons [19] with the five- 
dimensional domain- wall action for the fermions [20, 21]. 
Use of the Iwasaki gauge action in combination with 
domain-wall sea quarks allows for adequate tunneling 
between topological sectors [22], and in combination 
with domain-wall valence quarks reduces chiral symme- 
try breaking and the size of the residual quark mass as 
compared to the Wilson gauge action [23]. 

We compute the RHQ &-quark parameters on several 
ensembles with different light sea-quark masses; this al- 
lows us to study the sea-quark mass dependence, which 
we find to be statistically insignificant. We also deter- 
mine the parameters at two lattice spacings; we refer to 
the coarser ensembles with a w 0.11 fm as the "24 3 " en- 
sembles and the finer ensembles with a ~ 0.086 fm as the 
"32 3 " ensembles. Use of two lattice spacings allows us to 
take a naive continuum limit of physical quantities such 
as meson masses and splittings, although wc still include 
a conservative power-counting estimate of the residual 
C(|pa| 2 ) discretization errors from the RHQ action that 
may not be removed with this approach. Table I shows 
the parameters of the ensembles used for the RHQ pa- 
rameter tuning and bottomonium spectroscopy presented 



in this work. On the finer lattice spacings we double the 
statistics by performing two fermion inversions per gauge 
configuration with the origins of the quark sources sepa- 
rated by half of the temporal lattice extent. 

The ensembles listed in Table I have already been uti- 
lized to study the light pseudoscalar meson sector; we 
can therefore take advantage of many results from this 
earlier work. The amount of chiral symmetry breaking 
in the light-quark sector can be parameterized in terms 
of an additive shift to the bare domain-wall quark mass 
called the residual quark mass. At the values of M 5 = 1.8 
and L s — 16 used by RBC/UKQCD, the size of the resid- 
ual quark mass is quite small; am les — 0.003152(43) on 
the 24 3 ensembles and am TCS = 0.0006664(76) on the 32 3 
ensembles [10]. In order to compute the masses of B s 
and B* mesons for the tuning procedure we also need 
the value of the physical strange-quark mass on these en- 
sembles. This was already determined in Ref. [10]; am s — 
0.0348(11) on the 24 3 ensembles and am s = 0.0273(7) on 
the 32 3 ensembles. (In practice we use slightly different 
values of the strange-quark mass — am s = 0.0343 on the 
24 3 ensembles and am s = 0.0272 on the 32 3 ensembles — 
because we began this work before the light pseudoscalar 
meson analysis in Ref. [10] was finalized. These values, 
however, are within the stated statistical errors.) Finally, 
we must convert lattice meson masses into physical units 
for the tuning procedure and for comparison between 
predictions and experiment. The lattice scale was de- 
termined from the mass to be a -1 = 1.729(25) GeV on 
the 24 3 ensembles and a" 1 = 2.281(28) GeV on the 32 3 
ensembles [10]. These values are consistent with an inde- 
pendent determination of the 24 3 and 32 3 lattice spacings 
using the T(2S) - T(15) mass-splitting by Meinel [24]. 



B. Heavy-light meson correlator fits 

We extract the B s and B* meson energies from the 
exponential behavior of the following 2-point correlation 
functions: 

C Bs (t,t ;p) = J2 e%P ' S (oUy^)Op(0,h)) , (16) 
y 

C B: (t,t ;p) = l -Y,Y,^-y{0^{y,t)O v ^M)) , (17) 

i y 

where Op and Oy t are the pseudoscalar and vector 
heavy-strange meson interpolating operators, respec- 
tively: 



Op = 675s, Vi = tryis, 



(18) 



and the index "i" denotes the three spatial directions. We 
will explain the meaning of the tilde on some of the op- 
erators in Eqs. (16) and (17) later in this section. At suf- 
ficiently large times, excited-state contributions to these 
correlators will die away and the correlators will fall off 
as an exponential function of the meson ground-state en- 
ergy exp[— E(p)(t — to)]. We can therefore obtain the 
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TABLE I: Lattice simulation parameters used in our determination of the RHQ parameters for b-quarks and in our predictions 
for the bottomonium masses and mass-splittings. The columns list the lattice volume, approximate lattice spacing, light (mi) 
and strange (m^) sea-quark masses, unitary pion mass, and number of configurations and time sources analyzed. 
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ground-state energy from the following ratio of correla- 
tors: 



E c fr(p) = lim cosh 1 

t»t 



C(t,t ;p) + C(t + 2,t ;p) 



2C(t + l,t ;p) 



(19) 



which we refer to as the "effective energy" . In the above 
equation and throughout the remainder of this work, me- 
son masses and energies are given in lattice units (where 
the factor of "a" is implied) unless other units (e.g. GeV) 
are specified. 

We use the Chroma lattice QCD software system to 
compute the heavy and strange-quark propagators, as 
well as the 2-point correlation functions [25]. In order 
to minimize autocorrelations between data on nearby 
configurations, we translate the gauge field by a ran- 
domly chosen 4-dimensional vector before computing the 
strange-quark and 6-quark propagators. We generate the 
domain- wall light-quark propagators with a local (point) 
source; this allows them to be re-used for a future com- 
putation of B-meson decay constants and mixing ma- 
trix elements. In order to suppress excited-state con- 
tamination we generate the 6-quark propagators with a 
gauge- invariant Gaussian source for the spatial wavefunc- 
tion [26, 27]: 



b(x,t) = ^5(f,y;c7,JV)6(y,t) 



(20) 



where the smearing function S(x, y) depends upon the 
width a and the number of smearing iterations N: 



S($,y*;a,N)=[l + —V%g' , 



N 



(21) 
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As long as the parameters (a, N) satisfy the criteria N > 
3<r 2 /2, the source is spatially smooth and a good approx- 
imation to a Gaussian. For the free-field case (U = 1) 
with large N and small a, the root-mean-squared (rms) 
radius r rms ~ V3a/2 independent of N. Heavy-light me- 
son interpolating operators with a Gaussian-smeared b- 
quark are labeled with a tilde in Eqs. (16) and (17). We 



use a point sink, however, for both the strange and b- 
quark in the sink meson interpolating field because we 
find that this source-sink combination minimizes statis- 
tical errors in the correlators. 

Before beginning the iterative procedure to tune the 
RHQ parameters described in Sec. II B we compute the 
zero-momentum heavy-light meson pseudoscalar corre- 
lator [Eq. (16) with B s — > B{\ for several values of 
the Gaussian radius; these are given in Table II. Be- 
cause we expect both the light-quark and 6-quark mass- 
dependence of the optimal smearing choice to be mild, 
for each lattice spacing we analyze data on a single sea- 
quark ensemble and with a single light-quark mass and 
set of RHQ parameters {mod, cp,£}. For the smearing 
study on the 24 3 ensembles we use the preliminary results 
for the RHQ parameters in the chiral limit from Ref. [28] , 
{mod, cp, £} = {7.38, 3.89, 4.19}, which are similar to the 
earlier values presented at Lattice 2008 [11]. We analyze 
the unitary point on the am; = 0.005 ensemble. Fig- 
ure 2 shows the heavy-light pseudoscalar meson effective 
mass [E c g(p — 0)] for several choices of the Gaussian 
radius (including the limit of a point source) . The corre- 
lator generated with a 6-quark spatial wavefunction with 
a root-mean-squared (rms) radius of r rms = 0.777 fm 
clearly has the longest plateau with the earliest onset; 
we therefore choose to use this spatial wavefunction for 
the RHQ parameter tuning on the 24 3 ensembles. One 
might worry that the extremely long plateau in Fig. 2 is 
due to cancellations between excited states with positive 
and negative amplitudes, and does not correspond to the 
true ground-state mass. Figure 3 therefore shows a com- 
parison of the pseudoscalar meson effective mass in which 
the o-quark has a smeared source and point sink and one 
in which the 6-quark has both a smeared source and sink. 
The two effective masses agree within statistical errors, 
suggesting that we have obtained the true plateau. 

For the smearing study on the 32 3 ensemble we an- 
alyze the unitary point on the ami — 0.004 sea-quark 
ensemble. We use the RHQ parameters {moa,cp,C} = 
{3.70,3.60,2.20}, which are close to the preliminary re- 
sults on the ami = 0.004 ensemble in Ref. [12]. As in 
the case of the 24 3 ensembles, the Gaussian radius of 
''rms = 0.777 fm leads to the best plateau, so we use it 
for the RHQ parameter tuning procedure. This is con- 
sistent with expectations that the size of the £?-meson in 
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FIG. 2: Pseudoscalar meson effective mass for several choices for the Gaussian radius of the 6-quark in the heavy-light meson 
interpolating operator. Results are shown for the unitary point on the ami = 0.005 24 3 ensemble with RHQ parameters 
{m a,c P ,(} = {7.38,3.89,4.19}. 



TABLE II: Root-mean-squared radii and corresponding 
Chroma Gaussian smearing parameters [defined in Eq. (21)] 
considered here. The parameters shown in bold are used to 
obtain the RHQ parameters in the following subsection. 





a « 0.086 fm 


a « 0.11 fm 


r rms (fm) 


a 


N 


a N 


0.137 


1.39 


5 


1.83 5 


0.275 


2.78 


15 


3.6 25 


0.518 


5.24 


5 


6.92 80 


0.777 


7.86 


100 


10.36 170 


1.035 


10.48 


175 




1.047 






13.98 310 



physical units should be independent of the lattice spac- 
ing. 

We estimate the errors in the correlation functions and 
in the fitted meson energies using a single-elimination 
jackknife procedure. This allows us to propagate the 
statistical uncertainties including correlations between 
the parameters {mna,c p ,£} into subsequent steps of the 
RHQ parameter tuning procedure. We find no evidence 
of residual autocorrelations between subsequent trajecto- 
ries, as measured by comparing the errors between binned 
and un-binned data. We perform the v-squared mini- 
mization including the full covariance matrix, and choose 
fit ranges that yield acceptable correlated confidence lev- 



els (p- value 1 > 10%). 



Because the B s and B* meson energies are largely in- 
sensitive to the sea-quark mass, we expect the excited- 
state contamination to die off and the onset of the 
ground-state plateau to occur at around the same loca- 
tion on all sea-quark ensembles for a given lattice spacing. 
We therefore choose the same fit range for all sea-quark 
ensembles on a given lattice spacing. The requirement 
that we obtain a constant fit to the effective energy with 
a good correlated confidence level using the same fitting 
range for all ensembles helps to ensure that we obtain the 
true ground-state energy, and are not misled by "wiggles" 
in the plateau that are due to fluctuations in the gauge 
field, but are different on each ensemble. We do not, how- 
ever, expect excited-state contributions to be the same 
for all momenta, and, in fact, we observe an earlier onset 
for the plateau in the zero momentum effective energy 
than for the other momenta. Table III shows the fit- 
ting ranges used on the 24 3 and 32 3 ensembles. Figure 4 
shows the B s and B* meson effective energies for lattice 
momenta up to (ap) 2 — 3 on the ami = 0.005 24 3 en- 
semble. Effective energy plots for the other 24 3 and 32 3 
ensembles look similar. 



We adopt the PDG convention that the p- value is the probability 
of finding a \ 2 value greater than that obtained in the fit; hence a 
larger p- value denotes a stronger compatibility between the data 
and the fit hypothesis [29]. 
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FIG. 3: Pseudoscalar meson effective mass for the 6-quark Gaussian radius r rms = 0.777 fm. The full symbols correspond 
to correlators in which the 6-quark is generated with a Gaussian spatial wavefunction but has a point sink; the open points 
correspond to correlators in which the 6-quark has a Gaussian spatial wavefunction at both the source and sink. The effective 
masses agree, but the smeared-point data has smaller statistical errors. 



TABLE III: Time ranges used in plateau fits of the B a and 
Bg effective energies. We use different ranges for zero and 
nonzero momenta, but use the same range for all sea-quark 
masses at a given lattice spacing. 



fit range 

aw 0.11 fm [10,25] [10,25] 
aw 0.086 fm [11,21] [14,21] 



C. Determination of bottom-quark parameters 

We begin our iterative tuning procedure using the pre- 
liminary values for {moa,cp,£} determined in the pilot 
studies of Refs. [28] and [12]. We compute the B s and 
B* meson energies for seven sets of parameters surround- 
ing these values. We then determine the ratio of the 
rest mass to the kinetic mass for these seven parame- 
ter sets by fitting the nonzero momentum data for the 
B s meson to the energy-momentum dispersion relation, 
Eq. (5). Finally, wc determine the predicted values of the 
RHQ parameters from Eq. (15) using the experimentally- 
measured meson masses Mp B = 5.366 GeV and Mb* = 
5.415 GeV [29]. We find that the resulting values of 
{moa, cp,(} lie outside the "box" determined by the 
seven parameter sets. We therefore re-center the box 
around the newly-determined values and repeat the pro- 
cedure. We find that we need to iterate once or twice 
before the values of {mod, cp, C} settle down and remain 



TABLE IV: Final "box" of parameters used to obtain the 
tuned values of {moa, Cp, (} (see Fig. 1). In each column the 
first number is the central value of the parameter and the 
second number is the variation. 





moa 


cp 


c 


a w 0.11 fm 


8.40 ± 0.15 


5.80 ± 0.45 


3.20 ± 0.30 


a « 0.086 fm 


3.98 ± 0.10 


3.60 ± 0.30 


1.97 ± 0.15 



inside the box. Here we only show results for the final 
iteration, since plots for intermediate iterations look sim- 
ilar. The final sets of parameters used to obtain the tuned 
values of {moa, cp, C} on the 24 3 and 32 3 ensembles are 
given in Table IV. 

Figure 5 shows the energy-momentum dispersion rela- 
tion fit for both the B s and B* mesons on the ami = 
0.005 24 3 ensemble. Dispersion relation plots for the 
other sea-quark masses, RHQ parameter sets, and lat- 
tice spacing look similar. The slopes (M1/M2) of the 
B s and B* energy-momentum dispersion relations agree 
with unity (and hence with each other) within errors in 
the region of the parameter space near the tuned values 
of {moa, cp, £}. We choose to use the pseudoscalar me- 
son data, however, for the parameter tuning because it 
has smaller statistical errors. We perform a one parame- 
ter linear fit in which we fix the intercept to go through 
the measured value of the rest mass E(p = 0) and allow 
the slope to vary. We include data with lattice momenta 
through (ap) 2 — 3, and see no evidence for higher-order, 
e.g. 0([ap] 4 ), lattice discretization effects at these values 
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FIG. 4: Heavy-strange pseudoscalar meson (blue circles) and vector meson (red triangles) effective energies on the ami = 0.005 
24 3 ensemble with RHQ parameters {mod, cp, £} = {8.40, 5.80, 3.20}. From upper-left to lower-right the six plots show spatial 
momenta (ap) 2 = through (ap) 2 = 3. For each plot the shaded horizontal band shows the fit range used and the fit result 
with iackknife statistical errors. 



of the momenta. We account for correlations between 
data points by propagating the jackknife values of the 
energies from the 2-point fits described in the previous 
subsection. As a cross-check we compare the fit result 
with those of a two-parameter fit in which we allow both 
the slope and intercept to vary; we find that the results 
are consistent, and choose to use the one-parameter fit 
because it leads to smaller statistical errors in 

In order to reliably determine the RHQ param- 
eters via Eq. (15) we must be interpolating in a 
regime in which the bottom-strange meson observables 
{M Bs , AM Bs , Mf s /M 2 Bs } depend linearly upon the pa- 
rameters in the action {mod, cp,C}- We test this as- 
sumption and look for signs of curvature by comput- 
ing the observables for three different boxes of seven 
parameters with sizes ±a {moa ^ CpX} , ±2(T {moa ^ CpX} , and 
^ a {moa,c P ,C} (except for the parameter mpa on the 24 3 



ensemble for which the largest box is ±4cr moa ). We then 
determine the predicted values of the RHQ parameters 
for each of the three boxes; we find that the difference is 
negligible within statistical errors. 

Figure 6 shows the dependence of the spin-averaged 
mass, hyperfinc splitting, and rest mass over kinetic mass 
on the parameters moa, cp, and £j respectively, on the 
ami = 0.005 24 3 ensemble. We plot these dependencies 
because these are the parameters to which each observ- 
able is most sensitive. The bottom-strange observables 
{Mb b , AMp s , M^ s /M^" } depend linearly on the param- 
eters {moa, cp,£} throughout the range. The analogous 
plots for the other sea-quark ensembles look similar. 

Table V shows the nonperturbatively-tuned RHQ pa- 
rameters {moa, cp, C} obtained on the two 24 3 ensembles. 
We do not observe any statistically-significant sea-quark 
mass dependence. Hence, instead of extrapolating the 
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FIG. 5: B 3 (blue circles) and B* (red triangles) meson 
squared-energy difference versus spatial momentum-squared 
on the ami = 0.005 24 3 ensemble for the RHQ parameter val- 
ues {m a,cp,(} = {8.40,5.80,3.20}. The slope of the data 
gives the ratio of the meson rest mass over the kinetic mass 
(M1/M2). Data points shown with open symbol are not in- 
cluded in the fit. 



TABLE V: Tuned RHQ parameter values on the 24 3 ensem- 
bles determined using the parameter sets in Table IV. Be- 
cause we do not observe any statistically-significant sea-quark 
mass dependence, we obtain the final preferred values from an 
error-weighted average of the two sets of results. 



mo a 



cp 



c 



ami = 0.005 8.43(7) 5.7(2) 3.11(9) 
ami = 0.01 8.47(9) 5.8(2) 3.1(1) 
8.45(6) 5.8(1) 3.10(7) 



average 



RHQ parameters to the physical light-quark masses, we 
simply take an error-weighted average of the two values 
to obtain our final preferred results. Similarly, Table VI 
shows the nonperturbatively-tuned RHQ parameters on 
the three 32 3 ensembles and the corresponding weighted 
average. 



D. Comparison with perturbation theory 

It is useful to compare the nonperturbatively- 
determined values of the RHQ parameters with those 
computed in lattice perturbation theory. First, this pro- 
vides a consistency check of the nonperturbative tuning 
procedure. Second, this allows us to see how well the 
perturbative estimates are working in a case where we 
know the true nonperturbative value. Reasonable agree- 
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FIG. 6: Spin-averaged mass versus moa (upper plot), hy- 
perfine splitting versus cp (center plot), and rest mass over 
kinetic mass versus £ (lower plot) on the ami = 0.005 24 3 en- 
semble. The solid vertical lines with shaded gray error bands 
denote the tuned values of the RHQ parameters with jackknife 
statistical errors. For each quantity, the dashed line shows the 
dependence on moa, cp, or £ calculated from Eqs. (10)-(14). 
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TABLE VI: Tuned RHQ parameter values on the 32 3 en- 
sembles determined using the parameter sets in Table IV. Be- 
cause we do not observe any statistically-significant sea-quark 
mass dependence, we obtain the final preferred values from an 
error-weighted average of the three sets of results. 



moa cp ( 


ami = 
ami = 
ami = 


: 0.004 
: 0.006 
: 0.008 


4.07(6) 3.7(1) 1.86(8) 
3.97(5) 3.5(1) 1.94(6) 
3.95(6) 3.6(1) 1.99(8) 


averag 


;e 


3.99(3) 3.57(7) 1.93(4) 



merit between the two approaches bolsters confidence in 
our ability to rely on lattice perturbation theory in fu- 
ture situations where we do not have nonperturbative 
matching factors available. 

We calculate the RHQ parameters cp and £ at 1-loop 
in mean- field improved lattice perturbation theory [30]. 
The details of the perturbative calculation will be given 
in a separate publication [31]. The clover coefficient cp 
is obtained by matching the lattice quark-gluon vertex to 
the continuum counterpart in the on-shell limit. At in- 
termediate steps of the calculation infrared divergences 
are regulated with a nonzero gluon mass A; the final re- 
sults are obtained in the limit A — > 0. Similarly, the 
anisotropy parameter £ is obtained by requiring that the 
lattice heavy-quark dispersion relation, extracted from 
the momentum dependence of the pole in the heavy- 
quark propagator at one loop, agrees with the continuum. 
We implement the mean- field improvement in two ways. 
First we use the nonperturbative value of the fourth root 
of the plaquette, u = P 1 / 4 , to resum tadpole contribu- 
tions as in Ref. [30]. We also use the value of the spatial 
link field in Landau gauge to estimate uq. A compari- 
son of these two approaches is useful for ascertaining the 
systematic uncertainty due to the ambiguity in how to 
implement the tadpole resummation. The lattice per- 
turbation theory calculations of cp and £ also use the 
nonperturbatively-determined values of the bare-quark 
mass moa and the 2x1 rectangle R as inputs. The latter 
allows for a refined resummation of tadpole contributions 
in improved gauge actions [32]. 

Figure 7 compares results on the 24 3 ensembles in 
both unimproved and mean-field improved lattice per- 
turbation theory with the nonperturbatively-determined 
values. The results on the 32 3 ensembles look qualita- 
tively similar. The use of mean-field improved lattice 
perturbation theory brings the perturbative results into 
better agreement with the nonperturbative values. It 
also reduces the size of the one-loop corrections, thereby 
appearing to improve the convergence of the perturba- 
tive series, although one cannot be entirely sure that this 
trend persists to higher orders. In the case of cp, the 
unimproved one-loop corrections are very large (approx- 
imately a factor of 1.5) but are reduced to a more sensi- 
ble level by resumming tadpole contributions, whereas in 
the case of £ the unimproved one-loop corrections are 



already close to the naive power-counting estimate of 
a is? (V a 24 3 ) ~ 23% and the mean- field improved one- 
loop correctons are even smaller. 

We can use the results shown in Fig. 7 to estimate the 
uncertainties in the values of cp and £ calculated in lat- 
tice perturbation theory. We consider two approaches for 
obtaining the error. A naive power-counting estimate of 
the size of the neglected 2-loop corrections would lead to 
a predicted error of ct| ~ 5%. As mentioned earlier, how- 
ever, there is an ambiguity in how to estimate the tadpole 
factor Mo used in the resummation procedure. This is not 
strictly a measure of the size of higher-order corrections, 
but taking the difference between the values of cp and £ 
computed at one-loop using uq from the plaquette and 
from the spatial Landau link gives a larger estimate of the 
error in cp (~10-12.5%) than the naive power-counting 
approach. We therefore take this difference to be the er- 
ror in the perturbatively-calculated value of cp, but take 
Qj ~ 5% to be the error in the perturbatively-calculated 
value of £. For the central values we quote the average of 
the one-loop mean-field improved values expanded in the 
MS coupling at scale a -1 and computed with uq obtained 
from the plaquette and from the spatial Landau link. 

Our final perturbative estimates for cp and £ on 
the 24 3 and 32 3 ensembles are given in Table VII. 
They agree with the nonperturbatively-determined val- 
ues given in Table VIII in all cases. In order to pro- 
vide a fair comparison, we include an estimate of sys- 
tematic errors for both the perturbatively-calculated and 
nonperturbatively-computed values. The largest source 
of uncertainty in the lattice perturbation theory deter- 
minations is the error due to neglected terms in the 
coupling-constant expansion of 0(a 2 s ) and higher. In 
contrast, the largest source of uncertainty in the non- 
perturbative determinations of cp and £ is heavy-quark 
discretization errors from neglected operators in the ac- 
tion of 0(a 2 p 2 ) and higher (for moa the uncertainty in 
the lattice scale dominates). The good agreement be- 
tween lattice perturbation theory and the nonperturba- 
tive tuning procedure suggests that one-loop mean-field 
improved lattice perturbation theory is sufficiently reli- 
able that it can be used in situations where the nonper- 
turbative matching factors are not available, such as in 
our future computations of decay constants and mixing 
matrix elements. 



IV. BOTTOMONIUM MASS PREDICTIONS 

Given the determinations of the RHQ parameters de- 
scribed in the previous section, we can now make pre- 
dictions for other states involving 6-quarks, such as bot- 
tomonium masses and splittings. Comparison of the re- 
sults with experiment then provides a check of the rela- 
tivistic heavy quark framework and tuning methodology. 
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FIG. 7: Lattice perturbation theory calculations of cp (left plot) and £ (right plot) on the 24 3 ensembles [31]. From left to right, 
the perturbative calculations shown are (T) unimproved tree-level, (1) unimproved 1-loop, (Tp) mean-field improved tree-level 
using the plaquette to estimate the tadpole factor no, (Ipl), 1-loop mean-field improved value using the lattice coupling and 
the plaquette, (1pm) 1-loop mean-field improved value using the MS coupling and the plaquette, (Tu), mean-field improved 
tree-level using the spatial link in Landau gauge to estimate Uq, (Itjl), 1-loop mean-field improved value using the lattice 
coupling and the spatial link, and (Ium) 1-loop mean-field improved value using the MS coupling and the spatial Landau link. 
In each plot, the horizontal line indicates our choice of central value for cp or £ while the solid horizontal band denotes our 
estimate of the uncertainty with errors due to the truncation of the perturbative series and errors due to the uncertainty in 
mod added in quadrature. For comparison, the nonperturbatively-determined values are shown at the far right with statistical 
errors (solid inner error bar) and statistical and systematic errors added in quadrature (dashed outer error bar). 



TABLE VII: One-loop mean-field improved lattice perturba- 
tion theory predictions for the RHQ parameters cp and ( 
(right panel) [31]. The nonperturbative inputs used in the 
calculation - the bare heavy-quark mass mod, the plaquette 
P, the 2x1 rectangle R, and the spatial link in Landau gauge 
L - are given in the center panel. The errors in cp and £ are 
due to the truncation of lattice perturbation theory and the 
uncertainty in mod, respectively. The jackknife statistical er- 
rors in P, R, and L are negligible. 





nonperturbative 
inputs 
mod P R L 


perturbative 
estimates 

c 


a « 0.11 fm 
a w 0.086 fm 


8.45 0.588 0.344 0.844 
3.99 0.616 0.380 0.861 


4.8(6)(2) 3.2(2)(1) 
3.04(28) (7) 2.10(H)(5) 



A. Heavy-heavy meson correlator fits 

We extract the bottomonium meson masses from 
the following zero-momentum meson 2-point correlation 
functions: 



where is the bb meson interpolating operator for the 
state with spin structure T: 



£4 = &r& . 

DO 



(24) 



Table IX shows the interpolating operators used in 
the computation of the bottomonium 2-point functions. 
Again, the tilde over the interpolating operator in 
Eq. (23) denotes that the &-quark in the operator was 
generated with a Gaussian-smeared source. 

Plots of the effective energy, Eq. (19), for the bottomo- 
nium correlators show that excited-state contamination 
is significant for the choice of smearing that we used to 
obtain the RHQ parameters. In fact, on the 32 3 ensem- 
bles excited-state contamination appears to persist over 
the entire time range up to the temporal mid-point of 
the lattice, making a clean determination of the ground- 
state mass difficult. We therefore choose to use a dif- 
ferent smearing for the b quarks in the bottomonium 
correlators than for those in the bottom-strange corre- 
lators. We perform a similar smearing study to that 
described for bottom-strange states in Sec. IIIB. Fig- 
ure 8 shows the T (vector) and Xbo (scalar) meson effec- 
tive masses on the ami — 0.005 24 3 ensemble for several 
choices of the Gaussian radius and values of the RHQ 
parameters {m a,cp,(} = {8.40,5.80,3.20}. The corre- 
lator generated with a 6-quark spatial wavefunction with 
r [ms = 0.137 fm has the longest plateau with the earli- 
est onset; we therefore choose to use this spatial wave- 
function to compute the bottomonium masses and mass- 
splittings on the 24 3 ensembles. We perform an anal- 
ogous smearing test on the ami — 0.004 32 3 ensemble 
with RHQ parameters {m a,cp,(} — {3.70,3.60,2.20}. 
Again, we find that the Gaussian spatial wavefunction 
with r rms = 0.137 fm is best. Physically one expects a 
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TABLE VIII: Tuned values of the RHQ parameters on the 24 3 and 32 3 ensembles. The central values and statistical errors are 
from Tables V and VI. The systematic error estimates are obtained using the same approach as for the bottomonium masses 
and mass-splittings described in Sec. IV C. The errors listed in moo, cp, and C, are from left to right: statistics, heavy-quark 
discretization errors, the lattice scale uncertainty, and the uncertainty in the experimental measurement of the B s meson 
hyperfme splitting, respectively. Errors that were considered but were found to be negligible are not shown. For the scale 
uncertainty we quote smaller errors on the 32 3 ensembles because the lattice-spacing is determined more precisely than on the 
24 3 ensembles. 



moH cp £ 

;0.11fm 8.45(6)(13)(50)(7) 5.8(1)(4)(4)(2) 3.10(7)(11)(9)(0) 
■ 0.086 fm 3.99(3)(6)(18)(3) 3.57(7)(22)(19)(14) 1.93(4)(7)(3)(0) 



TABLE IX: Interpolating operators used to compute the 66 
2-point correlation functions. We average correlators over 
equivalent directions for the vector, axial-vector, and tensor 
states. 



meson 


operator 


Vb 


75 


T 


7i 


Xbo 


1 


Xbl 


7; 75 


hb 





bb meson to have a narrower spatial wavefunction than 
a bs meson, and this is consistent with our observa- 
tions. We find an optimal wavefunction that is approx- 
imately half as wide as the bottomonium rms radius 
r r ^ c s hardson = 0.224(23) fm computed from the Richard- 
son potential model [33] . 

Using the optimized 6-quark smearing, we then com- 
pute the bottomonium correlators, Eq. (23), on each en- 
semble for the final set of seven RHQ parameters used in 
the iterative tuning procedure. This enables us to propa- 
gate the statistical uncertainties in the RHQ parameters 
from the tuning procedure into our determinations of the 
bottomonium masses and mass-splittings. We determine 
the ground-state meson masses from constant fits to the 
effective mass. We observe similar excited-state contam- 
ination in the 775 and T states, so we choose a fit range 
that yields a good correlated confidence level for fits to 
both effective masses. Similarly, we use the same fit range 
for the Xbo, Xbi, an d /i& states. Finally, because we do 
not expect any significant sea-quark mass dependence, 
we use the same fit range for all sea-quark ensembles 
with the same lattice spacing. These constraints help to 
ensure that we are not fooled by false plateaus due to 
fluctuations in the gauge field, which will differ among 
uncorrelated ensembles. Table X gives the fit ranges to 
determine the various meson masses on the two lattice 
spacings. Figure 9 shows sample bottomonium effective 
masses and mass-splittings on the ami = 0.005 24 3 en- 
semble. Plots for other sea-quark ensembles (including 
at the finer lattice spacing) and other values of the RHQ 
parameters look similar. 



TABLE X: Time ranges used in plateau fits of the bottomo- 
nium effective masses. We use different ranges for the r/b and 
T states than for the x an d h states, but use the same range 
for all sea-quark masses at a given lattice spacing. 





fit range 


7]b & T 


XbO, Xbl, & hb 


aw 0.11 fm [15,30] 


[4,12] 


a « 0.086 fm [13,30] 


[7,20] 



B. Determination of bottomonium masses and 
fine-structure splittings 

We determine the predicted values of the bottomonium 
masses at the tuned RHQ parameters using equations 
similar to Eqs. (12)-(15): 
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(25) 



where the 1x3 matrix Jm and constant Am are de- 
termined from a finite difference approximation of the 
derivatives: 



Jm = 



M3 
2a 



M2 

moa 



M 5 - M 4 M 7 - M 6 
2er Cp 2dQ 



A M = Mi - J M x [m a, c P , (} 



(26) 
(27) 



and Mi is the bb meson mass measured on the i th parame- 
ter set listed in Eq. (11). (Note that the values of Mi, Jm, 
and Am are different for each bottomonium state.) For 
each jackknife set we use the values of the tuned RHQ pa- 
rameters {moa, cp, £} RH Q determined on that jackknife 
set, thereby preserving correlations between the three pa- 
rameters moa, cp, and £. Hence the jackknife statistical 
errors in the bb meson masses determined via Eq. (25) al- 
ready include the uncertainty due to the statistical errors 
in the tuned RHQ parameters. 

The use of Eqs. (26) and (27) requires that we are 
in a regime in which the bottomonium masses depend 
linearly on the RHQ parameters. We test this assump- 
tion and look for signs of curvature by computing the 
bottomonium masses for three different boxes of seven 
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FIG. 8: T (upper plot) and Xbo (lower plot) effective mass for several choices for the Gaussian radius of the fo-quark in 
the heavy-light meson interpolating operator. Results are shown for the ami = 0.005 24 3 ensemble with RHQ parameters 
{m a, c P , <} = {8.40, 5.80, 3.20}. 



parameters with sizes ±a {moa ^ CpX} , ±2a {moayCpX} , and 
i3<7/ mo (j,cp,C}- Figures 10 and 11 show the seven bot- 
tomonium masses and splittings — M rjb , Mr, My — M Vb , 
M X6 > M xm' M xti - M xbo> and M h h — versus m a, c P , 
and C on the ami = 0.005 24 3 ensemble; plots for the 
ami = 0.004 32 3 ensemble look similar. The statis- 
tical errors in the bb meson masses are approximately 
ten times smaller than those of the bottom-strange me- 
son masses, and we can resolve a nonlinear dependence 
of the bb meson masses on the RHQ parameters within 



statistical errors. This curvature is most pronounced in 
the hyperfme splitting Mr — M Vb , and the dependence 
is strongest upon the parameter £. The nonlinear de- 
pendence is mild, however, within the region of param- 
eter space defined by the inner-most box of parameters. 
Hence we expect that the use of Eq. (25) in this region 
will lead to only a small error in the bottomonium mass 
predictions. Nevertheless, we will include a systematic 
uncertainty in our predictions for the bottomonium me- 
son masses due to quadratic and higher-order corrections 
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FIG. 9: Bottomonium masses and mass-splittings on the ami = 0.005 24 3 ensemble with RHQ parameters {moa,cp,(} = 
{8.40, 5.80,3.20}. The meson states shown in each plot are specified in the legend. For each plot the shaded horizontal band 
shows the fit range used and the fit result with jackknife statistical errors. 
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to Eq. (25). 

Once we have the results for the bottomonium masses 
and mass-splittings at fixed sea-quark mass and lattice 
spacing, we must extrapolate to the physical light-quark 
masses and the continuum limit. Because the bb states 
contain no valence light quarks, we expect only a weak 
light-quark mass dependence and a correspondingly mild 
chiral extrapolation. In practice, as shown in Table XI, 
we do not observe any statistically significant dependence 
of the observables on the light sea-quark masses at either 
lattice spacing. We therefore compute the error- weighted 
average of each mass and mass-splitting over the different 
sea-quark ensembles at the two lattice spacings. 

Because the domain- wall fermion action is 0(a)- 
improved, the leading lattice discretization effects from 
the light-quark and gluon sector are proportional to 
a 2 . With the relativistic heavy-quark formalism, heavy- 
quark discretization errors depend on the lattice spacing 
as unknown functions of m^a [with coefficients of 0(1)] 
whose behavior is only known in the asymptotic limits 
of very large and very small moa; hence they do not 
have to scale as a 2 . As discussed in the following section, 
however, we estimate that gluon discretization errors in 
the bottomonium masses are larger than both light-quark 
and heavy-quark discretization errors, and consequently 
dominate the scaling behavior of the masses. We there- 
fore extrapolate the bottomonium masses to the contin- 
uum linearly in a 2 in order to remove gluon discretiza- 
tion errors. We estimate the remaining systematic un- 
certainty from heavy-quark discretization errors using 
power-counting, discussed below. Figure 12 shows the 
continuum extrapolation of the five bottomonium masses 
along with the experimentally-measured values for com- 
parison. 

In contrast, light-quark and gluon discretization errors 
largely cancel in the fine-structure splittings, so the scal- 
ing behavior is dominated by the heavy-quark discretiza- 
tion errors. With data at only two lattice spacings, how- 
ever, we cannot resolve quadratic or more complicated 
mo a dependence. We therefore choose not to extrapo- 
late the fine-structure splittings, and instead quote the 
results obtained on the finer 32 3 ensembles as our cen- 
tral values. Again, we estimate the residual systematic 
uncertainty from heavy-quark discretization errors using 
power-counting. 



C. Estimation of systematic errors 

We now discuss the sources of systematic uncertainty 
in the bottomonium masses and splittings. Table XII 
presents the total statistical and systematic error budget 
for each quantity. 



1 . Statistics 

We propagate the statistical errors through the en- 
tire multi-step analysis procedure via a single-elimination 
jackknife procedure. Hence the statistical errors include 
the uncertainty due to the statistical errors in the tuned 
RHQ parameters, including correlations between m^a, 
cp, and C- 

2. Heavy-quark discretization errors 

The RHQ action gives rise to nontrivial lattice-spacing 
dependence in physical quantities in the region moa ~ 1. 
Thus, instead of including additional functions of moa 
in the combined chiral-continuum extrapolation, we es- 
timate the size of discretization errors from the heavy- 
quark sector with power-counting. We follow the method 
outlined by Oktay and Kronfcld in Rcf. [34], in which 
they outline a general framework that applies to both 
heavy-heavy and heavy-light systems. 

We consider a nonrelativistic description of the heavy- 
quark action because both the lattice and continuum 
theories can be described by effective Lagrangians built 
from the same operators. Discretization errors arise due 
to mismatches between the short-distance coefficients of 
higher-dimension operators in the two theories. More 
precisely, for each operator Oi in the heavy-quark ef- 
fective Lagrangian, the associated discretization error is 
given by 

errorf Q = (C| at - C 4 cont ) (Of Q ) . (28) 

The "mismatch functions" fi = C\ at — C 4 cont are functions 
of the parameters of the lattice heavy-quark action. They 
have been calculated at tree-level for the anisotropic 
clover- improved Wilson action in Ref. [34]. The oper- 
ators d in Eq. (28) specify the 0(a 2 ) errors present 
in the heavy-quark action and their expectation values 
(Oi) depend on the physical quantity of interest. When 
the sizes of operators in the heavy-quark action are es- 
timated with power-counting appropriate to heavy-light 
meson systems, this framework leads to HQET. Similarly, 
when the sizes of operators in the nonrelativistic heavy- 
quark action are estimated with power-counting suitable 
for heavy-heavy meson systems, it leads to NRQCD. 

We consider two sources of heavy-quark discretization 
errors in the bottomonium system. The first is directly 
from operators that contribute to bottomonium masses 
and fine-structure splittings. The second is indirect con- 
tributions from discretization errors in the RHQ param- 
eters; these are due to heavy-quark discretization errors 
in the B s and B* energies used in the tuning procedure. 
We discuss each source briefly in turn and present the 
final error estimates here. Details are provided in the 
appendices A-C. 

To estimate the "direct" heavy-quark discretization er- 
rors, we compute the values of the mismatch functions for 
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FIG. 10: Bottomonium masses versus moo (upper plots), cp (center plots), and £ (lower plots) on the ami = 0.005 24 3 ensemble. 
The meson states shown in each plot are specified in the legend. The solid vertical lines with shaded gray error bands denote 
the tuned values of the RHQ parameters with jackknife statistical errors. For each quantity, the dashed line in the same color 
as the plotting symbol shows the dependence on the RHQ parameters calculated from Eqs. (25)-(27). 
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TABLE XI: Bottomonium masses and mass-splittings on the five sea-quark ensembles and averaged for each lattice spacing. 
For the masses, we extrapolate the results on the two lattice spacings to the continuum limit linearly in a 2 as described in the 
text. Errors shown are statistical only, but include the uncertainty due to the statistical errors on the tuned RHQ parameters. 







a 


« 0.11 fm 






a 0.086 fm 




continuum 


mass 


[MeV] 


ami = 0.005 


ami = 0.01 


average 


ami = 0.004 


ami = 0.006 


ami = 0.008 


average 




M Vb 




9328(14) 


9327(18) 


9328(11) 


9326(18) 


9341(15) 


9347(18) 


9338(10) 


9350(33) 


Mr 




9367(14) 


9367(17) 


9367(11) 


9379(16) 


9388(13) 


9395(16) 


9388(9) 


9410(30) 


M r - 


-M nb 


38.8(2.3) 


40.6(2.5) 


39.6(1.7) 


53.1(3.0) 


47.3(2.4) 


48.2(3.4) 


49.2(1.6) 




M Xb0 




9853(15) 


9848(18) 


9851(12) 


9816(19) 


9836(15) 


9837(20) 


9831(10) 


9808(35) 


M Xbl 




9884(15) 


9882(19) 


9883(12) 


9853(19) 


9873(15) 


9875(20) 


9868(10) 


9851(35) 


M Xbl 


~M Xb0 


31.2(1.8) 


33.5(2.0) 


32.3(1.3) 


37.8(2.7) 


36.6(2.2) 


38.8(2.6) 


37.5(1.4) 




M hb 




9895(16) 


9894(19) 


9895(12) 


9866(19) 


9884(16) 


9887(21) 


9879(10) 


9862(36) 



our lattice simulation parameters and estimate the sizes 
of the matrix elements of the higher-dimension operators 
Oi in Eq. (28) with power-counting appropriate to heavy- 
heavy meson systems. We use a -1 = 2.281 GeV [10], 
which is the lattice scale on our finer 32 3 ensembles, and 
m b = 4.2 GeV [29]. The RHQ parameters on the 32 3 lat- 
tices are given by {moa, cp,(} — {3.99,3.57,1.93}. We 
also need an estimate for the &-quark velocity v in the 
bb mesons. Following Ref. [2], we expect that the mass 
difference between the T(15) and T(25) states, which is 
roughly 500 MeV, should be of the same size as the aver- 
age kinetic energy, E ~ rribV 2 . Taking the quark mass to 
be half the meson mass gives an estimate for the 6-quark 
velocity squared of v 2 ~ 0.1. 

The numerical estimates of the relevant mismatch 
functions are given in Appendix A. Because the b quarks 
in the bb mesons are nonrelativistic, we estimate the size 
of operators using the "NRQCD" power-counting formu- 
lated in Ref. [3]: 

D ~ mt,v , gE ~ m b i r, gB ~ m b v , g ~ i> , (29) 

where the expansion parameter v is the spatial velocity 
of the b quarks. Thus, in NRQCD, an operator's numer- 
ical importance is determined by the order in the heavy- 
quark velocity v, rather than the dimension. Within the 
NRQCD power-counting framework, bb meson masses are 
approximately M ~ 2rrib, generic mass splittings such as 
Mr (25) — Mr (IS) are ~ mi,v 2 and fine-structure split- 
tings such as the hyperfine, spin-orbit, and tensor split- 
tings are ~ m^v 4 . 

In the RHQ approach we tune the coefficients of the 
dimension five operators in the Symanzik effective theory 
nonperturbatively; hence the leading heavy-quark dis- 
cretization errors come from operators of dimensions 6 
and 7 in the Symanzik effective theory (or alternatively 
the heavy-quark effective Lagrangian) that are omitted 
from the lattice action. The dominant errors in the bb me- 
son masses come from operators that are of 0(v 4 ) in the 
NRQCD power-counting. In Appendix B, we estimate 
the size of their contributions to bottomonium masses 
to be ~ 0.34%. Contributions from operators of 0{v A ) 
cancel in the fine-structure splittings, such that the dom- 
inant errors come from operators that are of 0(v e ). In 



Appendix B, we estimate the size of their contributions 
to hyperfine splittings to be ~ 32% and to x-state split- 
tings to be ~ 43%. The errors in the hyperfine splittings 
are smaller because they only come from operators con- 
taining the term a ■ B (and permutations thereof), where 
B is the chromomagnetic field. 

To estimate the "indirect" heavy-quark discretization 
errors from the bottom-strange mesons used in the RHQ 
tuning procedure, we use the same values of the mis- 
match functions but estimate the sizes of operator matrix 
elements with power-counting appropriate to heavy-light 
meson systems. We consider separately heavy-quark dis- 
cretization errors in the three input quantities: the spin- 
averaged rest mass Mb s , the hyperfine splitting AMb s , 
and the ratio of rest-to-kinetic masses M^ s /M^" ■ 

The 6-quarks in B hadrons typically carry a spatial mo- 
mentum \p\ w Aqcei, the scale of the strong interactions. 
Therefore we estimate the size of operators using HQET 
power-counting, which in the continuum is an expansion 
in \p\/mt,. The lattice introduces an additional scale, a. 
Following Ref [34], we therefore expand in powers of A, 
where A is either of the small parameters 

A ~ aA QCD , A QCD /m Q . (30) 

Within the HQET power-counting framework, bl meson 
masses are approximately M ~ to;, and hyperfine split- 
tings are ~ AQ CD /2m;,. 

As for the estimates above, we use the lattice-spacing 
and RHQ parameters on the 32 3 ensembles along with the 
experimentally-measured 6-quark mass. We also need an 
estimate for the 6-quark momentum Aqcei in the heavy- 
strange mesons. We choose Aqcd = 500 MeV because 
fits to moments of inclusive -B-decays using the heavy- 
quark expansion suggest that the typical QCD scale that 
enters heavy-light quantities tends to be larger than for 
light-light quantities [35]. 

The dominant errors in the B s and B* meson rest 
masses come from operators that are of 0{\ 2 ) in the 
HQET power-counting. In Appendix C, we estimate 

rC*) 

the size of their contributions to M 1 s to be ~ 0.05%. 
This is comparable to the size of the statistical errors in 
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FIG. 12: Continuum extrapolation of bottomonium masses and mass-splittings. Upper left plot: T (filled blue symbols) and 
rib (open red symbols) masses versus squared lattice spacing. Upper right plot: xt>i (open green symbols) and Xbo (filled pink 
symbols) masses versus squared lattice spacing. Lower plot: hb mass versus squared lattice spacing. On each plot the two 
lattice spacings a ~ 0.086 fm and a ~ 0.11 fm are indicated by vertical black dash-dotted lines. Data points at different 
light sea quark masses but the same lattice spacing are shown with an offset for clarity. The average values at each lattice 
spacing are given as shaded error bands in the same color as the symbols, and a linear extrapolation in a 2 of the averaged 
values leads to the continuum limit results denoted by circles. On the data points we show statistical errors only. On the 
continuum-extrapolated values we denote the statistical errors with solid error bars and the total statistical plus systematic 
errors with additional dashed error bars. For comparison we show the experimentally-measured values as stars. 



the effective masses computed in our numerical simula- 
tions (see the example fits in Fig. 4). As can be seen 
from Figs. 6, such a small variation in the spin-averaged 
mass leads to a statistically-negligible shift in the tuned 
value of mod (i.e. well within the vertical gray error 
band). Hence we neglect heavy-quark discretization ef- 
fects in M b s when determining the size of heavy-quark 
discretization errors in the tuned RHQ parameters. 

The dominant errors in the B s hyperfine splitting come 
from operators that are of C(A 3 ) in the HQET power- 
counting. In Appendix C, we estimate the size of their 
contributions to AMg s to be ~ 4.4%. This is approxi- 
mately twice as large as the statistical errors in the hy- 
perfine splittings computed in our numerical simulations. 
As can be seen from Figs. 6, a variation of this size leads 



to a statistically-significant shift in the tuned value of cp, 
so we must propagate it to an uncertainty in the tuned 
RHQ parameters. We estimate this error by varying the 
value of AMb s used in the RHQ parameter-tuning pro- 
cedure by ±4.4% and then re-computing the bottomo- 
nium masses and mass-splittings. For each mass or mass- 
splitting we take the largest variation observed on any of 
the sea-quark ensembles. We find that a ~ 4.4% error 
AMb s leads to a ~ 0.0-0.1% error in the bottomonium 
masses, a ~ 8.8% error in the hyperfine splitting, and a 
~ 6.2% error in the x-state splittings. 

Discretization errors in the B s kinetic meson mass arise 
from both the constituent quarks' kinetic energies and 
the binding energy. In Appendix C, we estimate their 
size to be ~ 2.6% following the method of Ref. [17]. This 
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is comparable to the size of the statistical errors in the B s 
meson kinetic masses computed in our numerical simula- 
tions (see the example fits in Fig. 5). As can be seen from 
Figs. 6, a variation of this size leads to a statistically- 
significant shift in the tuned value of £, so we must prop- 
agate it to an uncertainty in the tuned RHQ parameters. 
To estimate the resulting error we follow the same pro- 
cedure as described above for the discretization errors 
in the hyperfine splitting. We find that a ~ 2.6% error 
Aff'/Mf' leads to a - 0.1-0.2% error in the bottomo- 
nium masses, a ~ 3.6% error in the hyperfine splitting, 
and a ~ 1.0% error in the %-state splittings. 

To obtain the total heavy-quark discretization errors in 
the bottomonium masses and mass-splittings, we add the 
direct errors and the indirect errors in quadrature. The 
resulting estimates are given in Table XII. Numerically, 
the indirect errors due to discretization errors in the RHQ 
parameters turn out to be smaller than the direct errors 
for the 66-meson masses, and significantly smaller than 
the direct errors for the fine-structure splittings. 



3. Light-quark and gluon discretization errors 

We estimate the size of light-quark and gluon dis- 
cretization errors following the same approach as de- 
scribed for heavy-quark errors in the previous subsec- 
tion. In this case, the dimension 6 and higher-order 
light-quark and gluon operators in the Symanzik effec- 
tive Lagrangian have no counterpart in the continuum 
QCD Lagrangian. (There are no dimension 5 operators 
because both the light-quark and gluon actions are 0(a)- 
improved.) Thus the coefficients of the continuum opera- 
tors in the "mismatch functions" defined in Eq. (28) are 
£cont _ q F ur ther, the coefficients of the lattice oper- 
ators are not expected to be suppressed by any powers 
of the heavy-quark mass 1/toq. Thus we take them to 
be C] at = 1. The light-quark and gluon discretization er- 
rors are then given by expectation values of light-quark 
and gluon operators between heavy-heavy (QQ) meson 
states, i.e.: 

errorf Q ' 9 = (0^ Q ' 9 ) , (31) 

where we estimate their size using the NRQCD power- 
counting, Eq. (29). 

The largest discretization errors in bottomonium 
masses from the light-quark and gluon sector will arise 
from operators with only gluons. This is because any 
operators containing light-quark fields must extract light 
quarks from the sea, and their expectation values be- 
tween QQ meson states will be suppressed by at least a 2 . 
A typical dimension 6 gluon operator in the Symanzik ef- 
fective Lagrangian is 

Ogiue = trfF^D^] . (32) 



Within the NRQCD power-counting we expect its size to 
be 

(0 g i ue ) NR Q CD ~a 2 mV, (33) 

where two powers of mv come from the derivative oper- 
ators, and we estimate the size of F 2 to be the typical 
kinetic energy mv 2 . On the 24 3 (32 3 ) ensembles the cor- 
responding errors in the bottomonium masses are 

error gluo - a 2 m 3 v 4 /2m b = 3.0% (1.7%) , (34) 

which are several times larger than the estimated sub- 
percent contributions of heavy-quark discretization er- 
rors. Thus we conclude that, for bottomonium masses, 
the 0(a ) light-quark and gluon discretization errors will 
dominate the scaling behavior, and we can remove them 
by extrapolating to the continuum limit in a 2 . The sta- 
tistical errors in the continuum-limit values reflect the 
uncertainty on the slope in a 2 . 

Contributions from light-quark and gluon operators 
will largely cancel in the bottomonium fine-structure 
splittings, and we expect their contributions to these 
quantities to be negligible as compared to the heavy- 
quark discretization errors estimated previously. 

4- Input strange-quark mass 

We tune the parameters of the RHQ action from the 
bottom-strange system using the determination of the 
bare strange-quark mass on the two lattice-spacings from 
RBC/UKQCD's analysis of the light-pseudoscalar me- 
son sector in Ref. [10]. Hence the uncertainty in the 
bare strange-quark mass leads to a systematic error in 
the RHQ parameters, and consequently in the bottomo- 
nium masses and mass-splittings. We estimate this er- 
ror by varying the valence strange-quark mass in the B s 
and B* meson correlators used for the tuning procedure, 
Eqs. (16) and (17), and then re-computing the bottomo- 
nium masses and mass-splittings. 

Figure 13 shows the dependence of the meson masses 
and mass-splittings on the valence strange-quark mass 
used to tune the parameters of the RHQ action on the 
ami — 0.005 24 3 ensemble. The results at the three 
strange-quark mass values are consistent within statis- 
tical error, and analogous plots on the ami = 0.004 32 3 
ensemble look similar. Because the s» 1.2% uncertainty 
in m s leads to a 0.1% or less change in the bottomonium 
masses and a 0.3% or less change in the mass-splittings, 
we can safely neglect its effect from our error budget. 

5. Input scale uncertainty 

At first glance, the value of the lattice spacing in phys- 
ical units enters the computation of the bottomonium 
masses and mass-splittings in two ways. It first en- 
ters indirectly through the parameters of the RHQ ac- 
tion, which we tune by matching the values of the B s 
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FIG. 13: Bottomonium masses and mass-splittings versus the 
valence strange-quark mass in the bottom-strange meson cor- 
relators used to tune the parameters of the RHQ action. Re- 
sults are shown for the ami = 0.005 24 3 ensemble. The meson 
states shown in each plot are specified in the legend. For each 
quantity, the thicker line in the same color as the plotting 
symbol is an uncorrelated linear fit used to obtain the slope 
AAfa/Am,. The vertical solid line with gray error band de- 
notes the value of the physical strange-quark mass obtained 
in Ref. [10]. For each quantity, the two horizontal dashed 
lines show where the linear fit crosses the edges of the error 
band, thereby indicating the error due to the uncertainty in 
the strange-quark mass. 



and B* meson masses obtained on the lattice to the 
experimentally-measured values from the PDG [29]. It 
then enters directly when we convert the lattice values of 
the bottomonium masses and mass splittings into GeV 
in order to compare with experiment. In fact, however, 
the RHQ parameter tuning procedure allows us to avoid 
this second source of scale uncertainty. This is because 
our lattice calculation of the mass Mr* of a bb meson 
gives directly the dimensionless ratio M bb /MB S at the 
tuned values of the RHQ parameters {m a, cp, C} with- 
out further reference to the lattice scale. By construc- 
tion, at the tuned point the i? s -meson mass is fixed to 
the experimentally-measured value; hence we can pre- 
cisely obtain the bottomonium mass or mass-splitting in 
GeV by multiplying the ratio by M Bs = 5.3366 GeV [29]. 
We therefore need only consider the implicit dependence 
on the lattice spacing due to the RHQ parameters when 
estimating the scale uncertainty in the 66-meson masses. 

The absolute lattice scale (a -1 ) has a quoted statisti- 
cal error of - 1% (1.5%) on the 32 3 (24 3 ) lattices [10], 
where the errors on the two lattice spacings are highly 
correlated because they come from a single fit to data 
at both lattice spacings. We estimate the corresponding 
error in the bottomonium states by varying the lattice 
scale a -1 used in the RHQ parameter tuning procedure 
by plus and minus a statistical sigma on each sea-quark 
ensemble. For each bottomonium mass or mass-splitting, 
we then take the largest variation on any of the ensem- 
bles to be the error due to the uncertainty in the lattice 
scale. We find that the resulting uncertainty in the me- 
son masses is 0.2% or less, and in the mass-splittings is 
~l-3%; these errors are given in Table XII. 



6. Experimental inputs 



We tune the parameters of the RHQ action by using 
the experimental measurements of the spin-averaged B s 
meson mass and hyperfine splitting. The B s and B* me- 
son masses are both known to sub-percent precision [29] , 
so the experimental error in Mg s contributes a negligi- 
ble uncertainty to the tuned values of the RHQ param- 
eters. The experimental error in the hyperfine splitting 
AM Bs = 49.0(1.5) MeV [29], however, is - 3.1% and 
cannot be neglected. We estimate the error in the bot- 
tomonium masses and mass-splittings due to the exper- 
imental uncertainty in the B s meson hyperfine splitting 
by varying the value of AMb s used in the RHQ tuning 
procedure by plus and minus 1.5 MeV. For each bottomo- 
nium mass or mass-splitting, we then take the largest 
variation on any of the ensembles to be the correspond- 
ing error. We find that the resulting uncertainty in the 
meson masses is 0.1% or less, and in the mass-splittings 
is ~4-6%; these errors are given in Table XII. 
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7. Linear approximation 

We interpolate to the tuned values of the RHQ param- 
eters assuming a linear dependence upon {mod, cp,(}- 
Hence any deviation from linearity must be accounted 
for in the systematic error budget. In practice, as shown 
in Figs. 6, we do not see any statistically significant de- 
viation from linearity for the heavy-strange states over 
a wide range of RHQ parameters. Nor do we observe 
any statistically significant curvature for the x states or 
the hb (see the right-hand plots in Fig. 10). Thus the 
systematic uncertainty in the x states and the hb due 
to nonlinear dependence upon the RHQ parameters is 
negligible. We can resolve nonlinear dependence of T 
and rjb meson masses and the hyperfine splitting within 
the statistical errors in the measured effective masses, as 
shown in Figs. 10 and 11. The statistical errors in these 
data points, however, are almost two orders of magni- 
tude smaller than the statistical errors in the T and r\b 
meson masses and the hyperfine splitting interpolated to 
the tuned RHQ parameters given in Table XI; this is be- 
cause the interpolated values include the uncertainty due 
to the statistical errors in {mod, cp, (}. Hence we con- 
clude that the systematic error due to deviations from 
linearity is negligible for all bottomonium quantities con- 
sidered here. 



V. RESULTS AND CONCLUSIONS 

The relativistic heavy-quark formalism enables the de- 
scription of systems involving 6-quarks, such as B-mesons 
and bottomonium states, on currently available lattice 
spacings with lattice discretization errors from the heavy- 
quark sector of the same size as those from the light- 
quark sector. We have determined the 6-quark parame- 
ters for the RHQ action on the RBC/UKQCD 2+1 flavor 
domain- wall lattices with lattice spacings a w 0.11 fm 
and a ~ 0.08 fm. This is a continuation of and improve- 
ment upon the work of Li and Peng, who each presented 
preliminary results for B-mesons and bottomonium at 
conferences [11, 12]. 

In this work we tune the three parameters {moa, cp, C} 
using the bottom-strange system, where discretization er- 
rors are expected to be of 0([pa] 2 ) with \p\ Aqcd- 
We obtain the parameters nonperturbatively by impos- 
ing three simple conditions: that the masses of the B s 
and B* mesons agree with the experimental measure- 
ments, and that the B s meson on the lattice obey the 
continuum relativistic dispersion relation E 2 = p 2 + M 2 . 
We then test the reliability of the tuned parameters and 
the validity of the relativistic heavy-quark approach by 
making predictions for the masses and mass splittings of 
several bottomonium states. 

As shown in Fig. 14 and Table XIII, we obtain bot- 
tomonium masses with ^0.5-0.6% total uncertainties 
and mass-splittings with ^35-45% uncertainties, and 
find good agreement between our predicted values and 



experiment for all the quantities that we study. In fact, 
the preliminary work of Li successfully predicted the 
mass of the hb meson [11] before it was first observed 
by the Belle collaboration [36], thereby lending further 
credence to the relativistic heavy-quark formalism. We 
also find agreement with calculations of M Vh , the hyper- 
fine splitting M r - M m , and M hb using the NRQCD for- 
malism for the 6-quark [24, 37] and with a calculation of 
the hyperfine splitting using the Fermilab formalism [38] . 
Both the HPQCD and Fermilab/MILC works use the 
MILC collaboration's gauge configurations with 2+1 fla- 
vors of Asqtad- improved staggered sea quarks [39]; our 
study of bb meson spectroscopy using three flavors of dy- 
namical domain-wall light quarks provides a fully inde- 
pendent check of these results. Although the calculation 
by Meinel [24] uses the same RBC/UKQCD domain- wall 
+ Iwasaki configurations as in this paper, our result is 
still largely independent of his work because statistical 
errors (which are somewhat correlated between the two 
results) are not the primary source of uncertainty. 

Given the successful predictions of the bottomonium 
states, we now plan to use the nonperturbatively tuned 
parameters of the RHQ action to calculate _B-meson weak 
matrix elements of interest to flavor physics phenomenol- 
ogy. We are currently computing the leptonic decay con- 
stants fg and fg s and the neutral B°-B° mixing parame- 
ters [43] . These calculations are particularly timely given 
the observed approximately 3cr tension in the CKM uni- 
tarity triangle [44-47] which currently favors the presence 
of new physics in B^-mixing or B — > rv decays. Even- 
tually we would also like to use the RHQ framework to 
calculate more challenging quantities such as B — > n£v 
and B — > D^*Hv scmileptonic form factors, which are 
needed to extract the CKM matrix elements \V u b\ and 
\V c b\, respectively, from exclusive channels. Like the Fer- 
milab interpretation, our relativistic heavy-quark formal- 
ism applies to any value of the quark mass, and allows for 
a continuum limit. (This is in contrast to the NRQCD 
formalism, for which errors increase away from the in- 
finite heavy-quark limit.) Hence the same framework 
can be used for charm quarks, which are neither partic- 
ularly heavy compared to Aqcd nor light enough to be 
treated with a standard lattice light-quark formulation 
with 0(m c a) 2 errors that are well-controlled. Treatment 
of both b- and c-quarks within the same framework al- 
lows for further tests of the methodology. We therefore 
also plan to tune the parameters of the relativistic heavy- 
quark action for charm quarks, such that we can compute 
the leptonic decay constants fp> and fo s , as well as other 
weak matrix elements such as the short-distance contri- 
bution to D°-D° mixing. 

This work demonstrates the validity of the relativis- 
tic heavy quark action on bottom systems and opens 
a practical approach to obtain bottom and charm weak 
matrix elements with high precision given the computer 
resources currently available. Lattice QCD calculations 
of heavy-light weak matrix elements provide critical in- 
puts to the CKM unitarity triangle analysis. Hence de- 
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TABLE XII: Error budget for bottomonium masses and mass-splittings. The estimates of the size of each systematic uncertainty 
are given in the main text. Each error is given as a percentage, and we obtain the total systematic by adding the individual 
systematic uncertainties in quadrature. Errors that were considered but were found to be negligible (i.e. light-quark and gluon 
discretization errors, strange-quark mass uncertainty, and linear approximation) are not shown. 





M Vb 


Mr 


M r -M Vb 


M Xb0 


M Xbl 


M Xbl -M Xb0 


M hb 


statistics 


0.4 


0.3 


3.3 


0.4 


0.4 


3.7 


0.4 


heavy-quark discretization errors 


0.4 


0.3 


33.0 


0.4 


0.4 


43.6 


0.4 


input scale uncertainty 


0.2 


0.2 


3.2 


0.1 


0.1 


1.0 


0.1 


experimental inputs 


0.0 


0.1 


6.2 


0.0 


0.0 


4.3 


0.0 


total systematic 


0.4 


0.4 


33.7 


0.4 


0.4 


43.8 


0.4 



TABLE XIII: Comparison of predicted bottomonium masses and mass-splittings with experiment and, where possible, with 
other 2+1 flavor lattice calculations. The HPQCD and Meinel calculations use the NRQCD action for the 6-quarks [3], while 
the Fermilab/MILC calculation uses the Fermilab action [4]. For our results, the first error is statistical and the second is 
systematic; for the other results we add the errors in quadrature and quote the total. All results are given in MeV. 



this work Experiment HPQCD [40] Fermilab/MILC [38] Meinel [24] 

~M^ b 9350(33) (37) 9390.9(2.8) [29] 9390(9) 9400.0(7.7) 

Mr 9410(30) (38) 9460.30(26) [29] 

Mr-M Vb 49(02) (17) 69.3(2.8) [29] 70(9) 54.0(1^-|) 60.3(7.7) 

M Xb0 9808(35)(39) 9859.44(52) [29] 

M Xbl 9851(35)(39) 9892.78(40) [29] 

M Xbl -M Xb0 38(01) (16) 33.3(5) [41] 

M hb 9862(36)(39) 9899.1(1.1) [42] 9905(7) 9899.8(1.0) 



terminations with a variety of methods and indepen- 
dent sources of systematic uncertainty will be essential 
to definitively uncovering new physics in the flavor sec- 
tor. Use of the relativistic heavy-quark formalism for 
6-quarks on the RBC/UKQCD dynamical domain- wall 
lattices will provide phcnomcnologically- important, inde- 
pendent determinations of key heavy-light weak-matrix 
elements with comparable errors to other methods. 
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FIG. 14: Comparison of predicted bottomonium masses (left 
panel) and mass-splittings (right panel) with experiment. For 
the bottomonium masses we extrapolate the results on the two 
lattice spacings to the continuum linearly in a 2 , whereas for 
the fine-structure splittings we take the results on the finer 
32 3 ensembles as our central value. The solid error bars on the 
data points show the statistical errors. For our preferred re- 
sults, we also show the systematic errors added in quadrature 
as dashed error bars. 



Appendix A: Heavy-quark mismatch functions 

In this section we collect the forms of the mismatch 
functions used to estimate the size of heavy-quark dis- 
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cretization errors in heavy-heavy and heavy-light systems 
for the RHQ action. 

For each operator in the heavy-quark effective La- 
grangian, the "mismatch function" is defined as the dif- 
ference between the short-distance coefficients in the lat- 
tice and continuum theories. Hence the mismatch func- 
tions depend upon the parameters of the lattice action. 
The mismatch functions have been calculated at tree- 
level for the anisotropic clover-improved Wilson action 
in Ref. [34], but we present them here for completeness. 
Although Oktay and Kronfeld derive general expressions 
for ce 7^ cb and r s ^ 1 and include dimension 6 and 
higher order operators in the lattice action, here we show 
the mismatch functions specific to the RHQ action. We 
obtain these expressions from those in Ref. [34] by setting 
ce = cb = cpjC^ and r s = f, and setting the coefficients 
of the dimension 6 and higher-order operators to zero. 

There are five relevant tree-level mismatch functions 
that enter our estimates of heavy-quark discretization er- 
rors. The first is 



/ B (m a,c P ,C) 



1 



f 



where 



1 



2C 2 



8m|a 2 Srn^a 2 ' 



c 



m 2 a moa(2 + moa) I + mod' 

i e cc P 



4m E a 2 [m a(2 + m a)] 2 m a(2 + m a) 



(Al) 

(A2) 
.(A3) 



The function Je vanishes when the "chromoelectric 
mass" uie equals the 6-quark's kinetic mass m 2 . The 
second tree-level mismatch function is 



fw 4 (m a,c P ,() = \w A , 
6 



where 



W4 



2( 2 



mo£t(2 + mod) 4(1 + moa) 



(A4) 



(A5) 



The short-distance coefficient w 4 multiples the Lorentz- 
symmetry violating pj term in the lattice 6-quark's 
energy-momentum dispersion relation; hence the mis- 
match function f Wi vanishes when W4 = 0. The third 
tree-level mismatch function is 



where 



/m 4 (m a, cp,C) 



8C 4 



1 



1 



8m|a 3 8m 2 a 3 ' 



(A6) 



m 3 a 3 



4C 4 + 8C 3 (l + m a) 

[moa(2 + m^a)] 3 [moa(2 + moa)] 2 

c 2 



(1 + m a) 2 



(A7) 



The short-distance coefficient m l a3 multiplies the (p 2 ) 2 
term in the 6-quark's energy-momentum dispersion re- 
lation, so the mismatch function / OT4 vanishes when 



TABLE XIV: Tree-level mismatch functions for the 
nonperturbatively-tuned parameters of the RHQ action on 
the 24 3 and 32 3 ensembles. 



}e fw 



fn 



f u. 



fr, 



a w 0.11 fm 0.0640 0.0499 0.0353 0.0505 0.0934 
a « 0.086 fm 0.0864 0.0681 0.0521 0.0596 0.1359 



m,4 = m 2 . The fourth tree- level mismatch function is 



where 



f w > B {m a,c P ,() = 



cp 



B • 



1 + m a 



(A8) 



(A9) 



The coefficient w' B leads to a spin-dependent contribu- 
tion to the lattice quark-gluon vertex, so the mismatch 
function f w i vanishes when w' B = 0. The fifth tree-level 
mismatch function is 



fm B , (m a,c P ,() 



4mi( 



where 



m,p,a 3 



1 e-(c P 

m|a 3 (1 -I- moa) 2 



(A10) 



(All) 



The function fm B , vanishes when 7714 = m 2 (as above) 
and cp — C- 

To estimate the size of heavy-quark discretization 
errors in our numerical simulations, we evaluate the 
mismatch functions in Eqs. (Al), (A4), (A6), (A8), 
and (A10) at the tuned values of the RHQ parameters 
given in Tables V and VI. For the 24 3 ensembles we use 
{moa, cp,C} = {8.45,5.8,3.10} and for the 32 3 ensem- 
bles we use {m a,cp,C} = {3.99,3.57,1.93}. The re- 
sults are presented in Table XIV. Because the size of the 
heavy-quark discretization errors is sensitive to the nu- 
merical values of the tree-level mismatch functions, we 
have also tried evaluating Eqs. (Al), (A4), (A6), (A8), 
and (A10) at the tree-level values of the RHQ param- 
eters {moa, cp,C}. We find that the results are similar 
to those in Table XIV. We therefore conclude that the 
mismatch functions given in Table XIV reflect the typi- 
cal size of such coefficients for our simulations, and use 
them for estimating the heavy-quark discretization errors 
in the following appendices. 



Appendix B: Discretization errors in heavy-heavy 
meson masses and fine-structure splittings 

In this section we estimate the size of heavy-quark 
discretization errors in heavy-heavy mesons and fine- 
structure mass-splittings using the framework described 
in Sec. IV C 2. To estimate the numerical size of the 
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operator matrix elements, we use the NRQCD power- 
counting given in Eq. (29), and for the size of the coeffi- 
cients we use the mismatch functions on the 32 3 ensem- 
bles given in Table XIV. 

1. Masses 

Here we consider operators of 0(v 4 ), which pro- 
duce the dominant discretization errors in bottomonium 
masses. Oktay and Kronfeld enumerate all dimension 
6 and 7 bilinear operators in the heavy-quark effective 
Lagrangian consistent with symmetries in Table III of 
Ref. [34] . We do not need to consider contributions from 
dimension 8 bilinears because they will be of 0(v 6 ) or 
higher. 

a. 0(a 2 ) errors 

There are two dimension six bilinears that are of 0(v 4 ) 
in the NRQCD power-counting: 

h{~f ■ D, a ■ E}h , (BI) 

h~/ 4 (D ■ E - E ■ D)h . (B2) 

The expected size of these operators is 

(0f;) NRQCD „ a 2 m 3 v 4 _ (B3) 

At tree level the coefficients of these operators are both 
equal to Eq. (Al). We therefore estimate the contri- 
bution to the error from each of these operators to be 

errors = f E {0 E ) NRQCD /2m b ~ 0.15% , (B4) 

where we obtain the relative error in the bb meson masses 
by dividing by 2rrib, the size of the meson masses in the 
NRQCD power counting. 

b. 0(a 3 ) errors 

There are two dimension seven bilinears that are also 
of 0(v A ) in the NRQCD power-counting: 

hDfh , (B5) 
h(D 2 ) 2 h (B6) 

and the expected size of these operators is 

(0 4 ) NRQCD ~<z 3 m^ 4 . (B7) 

At tree-level the mismatch function for the first operator 
is given by f W4 , Eq. (A4) , so we estimate its contribution 
to the error in bb meson masses to be 

errors = f Wi (0 4 ) NRQCD /2m b ~ 0.21% . (B8) 

The tree-level mismatch function for the second operator 
is given by f mi , Eq. (A6), so we estimate its contribution 
to the error in bb meson masses to be 

error m4 = f m4 (0 4 ) NRQCD /2to 6 ~ 0.16% . (B9) 



c. Total error 

We obtain the total heavy-quark discretization error in 
the bb meson masses by adding the errors from the differ- 
ent operators in quadrature, including O e twice because 
there are two such operators: 

error totai = (2 x error! + error, 2 „ 4 + error^ 4 ) 1/2 

- 0.34%. (BIO) 

2. Hyperfine splittings 

Only spin-dependent operators containing the term 
a ■ B where B is the chromomagnetic field (and permu- 
tations thereof), contribute to hyperfine splittings such 
as the mass difference M T — M rjb [48, 49]. There are five 
dimension 7 bilinear operators of this form in the heavy- 
quark effective action at 0(v 6 ): 

Ei&HDliViBiyh, (Bll) 

h{D 2 ,iY,- B}h, (B12) 

V ,. ,///>:, /),/>',/;,//. (B13) 

Try • Z?i£ B-y Dh, (B14) 

JiDiiS ■ BD,h . (B15) 

Only the first two operators in Eqs. (Bll) and (B12) 
have nonzero matching coefficients at tree- level [34] . The 
matching coefficients of the remaining three operators in 
Eqs. (B13)-(B15) are zero at tree-level [34], and have not 
been computed to one-loop. Higher-dimension operators 
in the heavy-quark effective Lagrangian such as h{D 2 , a- 
(D x E — Ex D)}h also contribute to hyperfine splittings 
at 0(v 6 ), but the full set of dimension 8 heavy-heavy 
bilinears has not been worked-out in the literature. 

Given our incomplete knowledge of the 0(v e ) bilin- 
ear operators and corresponding mismatch functions, we 
use a more naive error estimation procedure for the bot- 
tomonium hyperfine splittings. The leading contribution 
to the hyperfine splittings is ~ mi> 4 , so contributions of 
0(v 6 ) are suppressed by by a factor of v 2 ~ 0.1. Hence 
we expect that neglected 0(v 6 ) operators lead to 10% 
errors in hyperfine splittings. We can check this estimate 
for the two cases in which the mismatch functions are 
known, as shown below. 

a. 0(a 3 ) errors 

The expected size of the operators in Eqs. (Bll) 
and (B12) is 

(a. B ) NRQCD ~ a 3 m 4 b v 6 . (B16) 

The tree-level mismatch function for the first operator is 
given by f w > , Eq. (A8), so we estimate its contribution 
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to the error to be 

error^ = (0 ff . B ) NRQCI 7m 6U 4 ~ 3.72% , (B17) 

where we obtain the relative error in bb meson hypcrfinc 
splittings by dividing by mjo 4 , the size of the hyperfine 
splittings in the NRQCD power counting. The tree-level 
mismatch function for the second operator is given by 
fm B n Eq. (A10), so we estimate its contribution to the 
error in bottomonium hyperfine splittings to be 

error mB , = f mgl (O a . G ) miQCD /m b v 4 ~ 8.48% . (B18) 

Both of these estimates are consistent with the naive 
power-counting expectation of 10% based on the order 
in the &-quark velocity v. 

b. Total error 

There are five dimension 7 and an unknown number of 
dimension 8 operators in the heavy-quark effective action 
that contribute to the hyperfine splittings at 0(v 6 ) in the 
NRQCD power-counting. If we assume that there are the 
same number of 0(v 6 ) operators at dimensions 7 and 8, 
we arrive at the estimate 

1/2 

error^ff F = f 10 x (u 2 ) 2 ) = 31.62% . (B19) 

3. X" s t a t e splittings 

The fine-structure splitting between \ mesons (M Xbl — 
M Xba ) is a linear combination of the spin-orbit and tensor 
splittings: 

^in-orbit = 1 {5Mxb2 _ 2M ^ o _ 3Mxfci) (B2Q) 

A!? 801 ' = \ (3M Xtl - M Xb2 - 2M Xb0 ) . (B21) 

Hence it receives contributions from both the spin- 
dependent operators containing a ■ B considered above 
(which lead to the tensor splitting [48]) and from spin- 
dependent operators containing D x E where E is the 
chromoelectric field (which lead to the spin-orbit split- 
ting [49]). 

a. 0(v 4 ) errors 

There is one relevant bilinear at dimension 6 which is 
of 0(v A ) in the NRQCD power-counting: 

h{j ■ D, a ■ E}h . (B22) 

We estimate the size of its contribution to the error in 
the x-state splittings to be 

errors = / B (0 B ) NRQCD /m b w 4 ~ 29.30% . (B23) 



Note that the contribution of this operator to the x _ state 
splittings is not as large as the order in the 6-quark veloc- 
ity v would suggest because of the small numerical size 
of/ E . 



b. 0(v 6 ) errors 

The same 0(v 6 ) operators that contribute to the hy- 
pcrfinc splittings also contribute to the splitting between 
the x states. We therefore estimate their contributions 
to be the same size as for the hyperfine splittings: 

errors = 31.62% . (B24) 



c. Total error 

We obtain the total heavy-quark discretization error 
in the \ state splittings by adding the 0(v ) and 0(v 6 ) 
errors in quadrature, yielding 

error^f = (error 2 4 + error 2 6 ) 1/2 = 43.11% . (B25) 



Appendix C: Discretization errors in heavy-strange 
meson masses and hyperfine splitting 

In this section we estimate the size of heavy-quark dis- 
cretization errors in the heavy-strange meson quantities 
- the spin-averaged mass, hyperfine splitting, and ratio of 
rest-to-kinetic masses - used in the RHQ parameter tun- 
ing procedure. Again we use the framework described in 
Sec. IV C 2. To estimate the numerical size of the opera- 
tors, we use the HQET power-counting given in Eq. (30), 
and for the size of the coefficients we use the mismatch 
functions on the 32 3 ensembles given in Table XIV. 



1. Rest mass 

Because we tune the coefficients of the dimension 5 op- 
erators in the RHQ action nonperturbatively, the leading 
discretization errors come from operators of dimension 6 
and higher in the effective theory. There are two dimen- 
sion 6 bilinears of 0(X 2 ) in the HQET power-counting: 

h{j D,a E}h, (CI) 
hji{D ■ E - E ■ D)h . (C2) 

The estimated size of these operators is 

(O e )^ ~ « 2 A^ CD . (C3) 

We do not consider operators of dimension 7 and higher 
because they are all at least of C(A 3 ). At tree-level the 
coefficients of the operators in Eqs. (CI) and (C2) are 
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both given by Eq. (Al), so we estimate their contribu- 
tions to the error in the spin-averaged B s meson rest 
mass to be 



errors = f E (0 E 



HQET 



/Mi 



0.04% . 



(C4) 



By construction, we tune the RHQ parameters such that 
the spin- averaged rest mass equals the experimental value 
of \ (M Bg + 3M^ s ) , so we obtain the relative error in Mi 
by dividing by M B „ = 5.4028 GeV. We obtain the total 
heavy-quark discretization error in the spin-averaged B s 
meson rest mass by adding the contributions from the 
two operators in quadrature, which yields: 



Ml,B 

error totai 



or ~ 3 MeV. 



(2 x error|; 



1/2 



0.05% , 



(C5) 



2. Kinetic mass 

Discretization errors in the kinetic meson mass Mi 
arise from both the constituent quarks' kinetic energies 
and from the binding energy. The Appendix of Ref. [17] 
provides a semi-quantitative estimate of the discretiza- 
tion error in Mi (see also Ref. [50]). Although this es- 
timate is made assuming that both quarks in the meson 
are nonrelativistic, the result is interpreted a 'posteriori 
under the assumption that the strange quark is light and 
relativistic. We follow the same approach here. 

The tree-level discretization error in M 2 through 0{v 4 ) 
in the nonrelativistic expansion is given by [17] 



5M 2 = 



3m2 2 



- 1 



41(74(7712 a) 



(C6) 



where this result applies to 5-wave states. Note that the 
5M 2 is zero if the masses 777,4 = 7772 and the Lorentz- 
symmetry violating coefficient 7/J4 = 0. To estimate the 
numerical size of the discretization error in M 2 we replace 
(p 2 ) with Aq CD following the HQET power-counting pre- 
scription and use the expressions for 7772, 7774, and W4 
given in Eqs. (A2), (A5), and (A7). By construction, 
we tune the RHQ parameters such that the kinetic me- 
son mass equals the experimental value of the B s meson 
mass, so we obtain the relative error in M% by dividing 
by M Bs = 5.366 GeV. We obtain 



error 



M 2 ,B. 

total 



= 2.59% 



(C7) 



or ~139 MeV. 



3. Hyperfine splitting 



from the dimension 5 operator hiT, ■ Bh and is of 0(A) 
in the HQET power-counting. Because we tune the co- 
efficient of this operator nonperturbatively, there are no 
associated discretization errors. Thus we consider dis- 
cretization errors from operators of 0(A 2 , A 3 ). There are 
five dimension 7 bilinear operators of the type a ■ B in the 
heavy-quark effective action at 0(A 3 ); these are given in 
Eqs. (B11)-(B15). Operators of dimension 8 and higher 
in the heavy-quark effective Lagrangian are all of 0(A 4 ) 
or higher in the HQET power-counting. 



The expected size of the operators in Eqs. (Bll) 
and (B12) is 



<a. S ) HQET ~a 3 A 4 QCD . 



(C8) 



By construction, we tune the RHQ parameters such 
that we reproduce the experimental value of the bottom- 
strange hyperfine splitting M B -Mb,- Hence we divide 
the contributions of these operators by AM Bs — 49 MeV 
to obtain the relative error in the B s hyperfine splitting. 
The tree-level mismatch functions for the two operators 
are f w > [Eq. (A8)] and f mB , [Eq. (A10)], so we estimate 
their contribution to the error in the bottom-strange hy- 
perfine splitting to be 



error,. 



= /^(a. s > HQET /AM Bs 

~ 0.64%, (C9) 

= f mB ,{O a . B )^/AM Bs 

~ 1.46%. (C10) 



b. 0(a s a 3 ) errors 

The expected size of the operators in Eqs. (B13)-(B15) 
is also 



(O 



oB 



HQET 



x 3 A 4 



QCD 



(Gil) 



The mismatch functions of these operators, however, van- 
ish at tree-level [34]. Because they have not been com- 
puted to one-loop, we simply estimate their size to be 
a^ s ( 1/0323) = 0.22. Under this assumption, we esti- 
mate that the contribution of each of these operators to 
the bottom-strange hyperfine splitting is 



error as = a s (O a . B 



HQET 



/AM Bs ~ 2.36% . (C12) 



The bottom-strange hyperfine splitting receives con- 
tributions from spin-dependent operators containing the 
term a • B where B is the chromomagnetic field (and per- 
mutations thereof) [48, 49]. The leading contribution is 



This estimate is likely conservative, given that we would 
naively expect O{a s o?) errors to be smaller than O(o?) 
errors, due to the fact that we have not considered any 
possible suppression from the 1-loop mismatch functions. 
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Total error 



or - 2 MeV. 



We obtain the total heavy-quark discretization error 
in the bottom-strange hyperfine splitting by adding the 
errors from the different operators in quadrature, includ- 
ing error Qfl three times because there are three 1-loop 
operators: 



error 



AM b 

total 



error?,,, + errorfL 



4.40% , 



3 x error„ 



1/2 



(C13) 
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